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1 Introduction 

In this paper we investigate multiscale symmetric collocation approximation 
with Wendland compactly supported radial basis functions (RBFs) to solve 
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the Stokes problem 



vAu + Vp = f in Q, 
V ■ u = in Q, 
u = g on dfl, 



(1) 
(2) 
(3) 



where the region fi C R d , the viscosity v, f : Q — > R d and g : Q — > R d are 
given and we seek an approximate solution to the velocity u : fl — > R d and 
the pressure p : Q — > R. 

Radial basis functions (RBFs) have been increasingly important in the area 
of approximation theory. For solving partial differential equations (PDEs), 
RBFs with meshless collocation for PDEs have been investigated in [5JH] and 
for the Stokes problem in [22]. Matrix- valued, positive definite kernels have 
been studied in [TMMT1 IT21[T3"1I22] . Two excellent recent books covering practi- 
cal and theoretical issues related to RBFs are [4] and [21] . We recall that a func- 
tion if - : R d — > R is said to be radial if there exists a function ip : [0, oo) —> R 
such that tf'(x) = -0 ( 1 1 x| 1 2 ) for all x € R d , where |j • H2 denoting the usual Eu- 
clidean norm in R d . Then with a scaling factor S > 0, we can define a scaled 



A practical issue that arises is that of which scale to use for the radial basis 
functions. A small scale will lead to a sparse and consequently well-conditioned 
linear system, but at the price of the approximation power. Conversely, a 
large scale will have better approximation power but at the price of a poorly- 
conditioned linear system. 

Many examples may naturally exhibit multiple scales, for example, mod- 
elling fluid dynamics, where many different scales may be required. Of course, 
this comes at the price of having to select which scales to use in which regions 
but this is not the topic of this paper. 

The multiscale algorithm investigated in this paper is constructed over 
multiple levels, in which the residual of the current stage is the target function 
for the next stage, and in each stage, RBFs with smaller support and with 
more closely spaced centres will be used as basis functions. 

Such a multiscale algorithm for interpolation was first proposed in [5] 
and P2] but without any theoretical grounding. Theoretical convergence was 
proven in the case of the data points being located on a sphere [TO] and then 
extended to interpolation and approximation on bounded domains |23| . 

In [TT] we can find an analysis of multiscale algorithms for RBF collocation 
of elliptic PDEs on the sphere. 

The extension to considering the Stokes problem on a bounded domain 
changes the analysis significantly as matrix-valued kernels need to be con- 
sidered with divergence-free approximation spaces. The scaled kernel, conver- 
gence and stability analysis are significant new contributions. 
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Our approach differs significantly from the literature, where for example 
in [Tl] the authors consider the incompressible Navier- Stokes equations in its 
weak formulation, and then decompose the velocity into coarse and fine scales. 

In the next section we provide necessary background material regarding 
point sets and function spaces. Scction|3]describes our (single scale) symmetric 
collocation approximation and then Section 2] extends this to a multiscalc 
algorithm and provides proofs of convergence. Scction[S]providcs an analysis of 
the stability of the approximations. Section [6] provides numerical experiments 
to test the theoretical results. 



2 Preliminaries 

In this paper, we will use (scaled) compactly supported radial basis functions 
to construct multiscalc approximate solutions to the Stokes problem, that is, 
we form the solution over multiple levels. We will work with a given domain 
J? C R d . A kernel $ : Q x Q ->• K is also given. 

At each level, we will have a finite point set X C fl. We will define the 
mesh norm as 

hx,n ■= sup min ||x — Xj|| 2 , 
and the separation distance as 

Qx ■= imin|| Xj -x fc j| 2 , 

Z 3^k 

which arc measures of the uniformity of the points in X. Then for example, at 
each level i, we denote the mesh norm by hi. The selection of point sets with 
mesh norms decreasing in a specific way will form one of the requirements for 
convergence of our algorithms and the separation distance will be used for the 
stability analysis. 

We define the Sobolev spaces in the usual way. For a given domain, fl G R d , 
k € No, and 1 < p < oo, the Sobolev spaces Wp(fi) consist of all u with weak 
derivatives D a u £ L p (f2), \a\ < k. The semi-norms and norms are defined as 




D a u\\l p(a) I and \\u\\ w u {n) = I £ P^H^ 
/ \M<fc 



For p = oo, these definitions become 

\ u \w^{n) = sup \\ Dau \\L^(n) and ||u|| M /^ (f2) = sup \\D a u\\ Lao[n) . 

\ot\—k 
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Let 1 < p < oo, i G No, and < s < 1. Then we can define the fractional 
Sobolcv spaces W k+s (S7) as all u for which the two norms 



\ u \wp +s {n) 
ll u llw*+*(n) 

are finite. For the case p = 2, we write W$(f2) = H k (Sl) and L 2 (i7) = W 2 °(J2). 

The functions that we will be concerned with are defined on a bounded 
domain fl with a Lipschitz boundary. As a result, there is an extension operator 
for functions defined in Sobolev spaces which is presented in the following 
lemma jl] Theorem 1.4.5]. For further details, we refer the reader to [19]. 

Lemma 2.1 Suppose Q C M. d has a Lipschitz boundary. Then there is an ex- 
tension mapping E$ '■ H T (S7) — > H T (M. d ) defined for all non-negative integers 
r satisfying Egv\n = v for all v € H T {fl) and 

\\Esv\\ht(r*) < C\\v\\nT(Qy 

C will denote a generic constant. Since we also have ||w|| jj T (r2) < ||^sv||ff T (K £i )7 
this means that when we need to consider the H T (S7) norms of the errors at 
each level, we can carry out our error analysis in the H T (R d )-norm. This is 
advantageous, since we then have for g e H T (R d ) 

llfllli^CR*) = / l5M| 2 (l + IMl2) T du;, (4) 




1/p 



dx dy 



upon defining the Fourier transform as 

g(u) = (2n)- d/2 J ' g(x)e- ixT "dx. 



At each level, we will also require a kernel W : SI x SI — >• R. We will use 
the Wcndland compactly supported radial basis functions [21] with a (level- 
specific) scaling parameter 5 > 0. We recall that for a given spatial dimension 
d and smoothness parameter k £ N, the Wendland functions are defined as 



4>tA r ) 



i 



r{k)2 k - 1 
o 



s (1 - sf{s 2 - r 2 ) k - L ds for < r < 1, 
for r > 1, 



(5) 



with I 



LlJ 



1. This choice of I ensures the Wendland functions are 



positive definite. It is the reproducing kernel of a Hilbert space which is norm 
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equivalent to the Sobolev space H~5~ +k (R d ). For the Wendland basis func- 
tions, there exist two constants < c\ < c%, which depend on d and k, such 
that their Fourier transforms satisfy [3T] 

^(l + Miy^^^ <$ e , k (^)<c 2 (l + M 2 2 y^- k -\ u;eR d . (6) 
With a given kernel & and scaling factor 6 > 0, we define the scaled kernel 

as 

!P*(x) := <T<V ( M^V (7) 



Appropriate selection of the scaling parameters will also prove to be one 
of the important ingredients for convergence of our multiscale algorithm. 
We will need norm equivalence as stated in the following lemma from [2]. 

Lemma 2.2 For every 8 € (0, 8 a ) and for all g £ H T (M. d ), there exist con- 
stants < C3 < C4 such that 

C3\\g\\& s < h\\H-(R") < Ci5^ T \\g\\<p s . 

As we will be working with vectors, in particular for u, we will need to 
define vector-valued Sobolev spaces in the usual way as 

H T {Q) := H T {Q) x ... x H T (f2), 

with norm 

/ d \ 1/2 

l|f|| H ^):= fEll/iH^wj • ( 8 ) 

Now we define divergence- free approximation spaces in fl and in M. d . With the 
divergence of u : Q — > M. d defined as 

d 

V ■ u := ^ fifUj , 
i=i 

we define 

H T (i7; div) := {u g W{Q) : V ■ u = 0} , 

and 

H r (M d ; div) := J f € H r (M d ; div) : / ^TItP ( X + IMli) T+1 < 00 1 , 

1Mb 



with norm 



IflllH^ - ( 2 -)" d/2 / 4^#(i + HI^) r+1 d- 
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We note that H T (]R d ; div) is a subspacc of H T (K d ; div). We will also need that 
for 17 C R d being a simply connected domain with C^' 1 boundary for d = 2, 3 
and with t > 0, there exists a continuous operator 

E div : H T (17; div) -> H T (R d ; div), 

such that E div u|l7 = u for all u e H r (17;div) [22, Proposition 3.8]. This 
operator is defined as 

E div u := V x E s Tu 7 (9) 

where Eg is the extension operator defined in Lemma l2.1l and T is a bounded 
operator T : H r (17; div) ->• H r+1 (17) with r = k + 6 for fc e N and 6» £ [0, 1]. 

To measure the pressure, which is determined only up to a constant, we 
will use the norm 

||p||jjt(«)/h := inf |b + c||jiT(fl). 

cGR 

Finally, we will need the following "sampling" inequalities, which are valid 
for both scalars and vectors [T()UTTll2"2"] . 

Theorem 2.3 Let Q C R d be a bounded domain with Lipschitz boundary. Let 
t > d/2. Let X C Q be a discrete set having mesh norm h sufficiently small. 
For each w € H r '(17) with w\X = we have for < a < r that 

\\w\\ HH a) <Ch T -"\\w\\ HT{n) . (10) 

Theorem 2.4 Let t = k + s > d/2. Let 17 C M. d be a bounded domain having 
C k ' s smooth boundary. Let X C 917 be a discrete set with h sufficiently small. 
Then there is a positive constant C such that for all w € iJ T (17) with w\X = 
we have for < a < t — 1/2 that 

\M\h*wij) < Ch T - 1 '*-' r \\w\\ B rw (11) 

We also define a matrix- valued function $ : K d — > R" x ™ as being positive 
definite if it is even, so <&(— x) = 4?(x), symmetric, so ^(x) = $(x) T , and 
satisfies 

n 

^2 a j®( x j - Xfc)«fe > o, 

for all pairwise distinct Xj € ]R d and all e M" such that not all arc 
vanishing. 
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3 Symmetric collocation approximation 

We will first consider a single-scale approximant to the combined velocity and 
pressure vector v := (u,p) : R d -> R d+1 , following [TSIEIEI]. Then ©-© 
become 

d 

(Lv) i := djjVj + d l v d +i = U in fl, (12) 

3=1 

d 

J2 d 3 V J = in Q > ( l3 ) 

3=1 

Vi = g t on df2, (14) 

where 1 < i < d. We seek a meshfrce, kernel-based collocation method with an 
analytically divergence-free approximation space. We use the notation Vv+i 
and ipr-i to denote the functions to be used in our matrix-valued kernel. We 
will mainly be interested in the case where both and Vv-i arc Wcndland 
functions which, for a given spatial dimension d, have native space norms 
equivalent to the Sobolev spaces H T+1 (M. d ) and H T ~ 1 (M. d ) respectively. Their 
Fourier transforms satisfy 

ci, T+ i(l + |M||)- T - 1 < Vv+i(M| 2 ) < c 2)T+1 (l + llc^Hl)--- 1 , (15) 

and 

c w (l + ll^lll)-^ 1 < C(||o,|| 2 ) < c 2)T _ 1 (l + ||u,||l)- T+1 , (16) 

and we define C± := min(ci i7 -_|-i, Ci iT _i) and C 2 := max(c 2)T +i, c 2)T _i). Then 
we define the matrix- valued kernel 

* := V; ° J : R d -> R (.d+l)x(d+i) ) (17) 

where tf'-r+i := (-A1 + VV T )f/v+i with I denoting the identity matrix. We 
note that \P r +i is also positive definite [T5] and hence due to the tensor product 
construction of 3?, it is positive definite as well. This choice for \P T +i is known 
to lead to divergence- free interpolants [15] . We also note that 

C+7(w) = (|M|2l-^u; T )V^l(u;). (18) 

We will consider the case where the collocation points are the same as the 
RBF centres. We denote the interior centres by X\ := {xi, . . . ,xn} and the 
boundary centres by X 2 ■= {xn+i , ■ • ■ , xjj} and their union by X — X\ U X 2 , 
with mesh norms h\ and hi respectively. Since (|13|) is automatically satisfied, 
this means that our approximant and collocation conditions will consist of dN 
terms from (TT2"]) and d(M — N) terms from (TT4"]) . Then with L y denoting the 
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operator L acting as a function of the second argument, applied to rows of 
our approximant takes the form 

d N d M 

Sxv(x)^^ QlJ (L¥(x-x J )) i +^ ]T a,,,-* (x -*,-)«, (19) 

i=l j=l i=l j=N+l 



where the notation means column i of the matrix #. The coefficients 
1 < i < d, 1 < j < M arc determined by the collocation conditions 



(LSxv( Xj )). = /fa), 1 < i < d, j = 1, . . . , N, (20) 
(S x v(x J -)) i = 9i(xj), l<i<d, j = N+l,...,M. (21) 

From [6"ll22j. we know that if i/v+i? Vv-i are positive definite and if ^ T +i € 
W?(R d ) n C 2 (R d ), then the native space of the kernel * given by dTTJ) is 



A/"*(R d ) = A/* T+1 (R d ) x M^_ 1 (R d ), 



with norm 



f llA^.(» d ) 



IfuIlL 



(2tt 



-d/2 



|f u («)||i , |/ p (u,)| 2 



|w||iVr+l(w) ^r-l(w) 



do;, 



(22) 



where f = (f u ,/ P ) T with f u 



and / p 



We recall that the 



generalised interpolant satisfies Chapter 16] 

||v - SxEvH^Rd) < Hvlljv^^d). 

With (|15[) and (|16|) . upon defining the extension operator for the velocity- 
pressure vector v as 

Ev := (E div u,£ S p) , (23) 



where E$ is the classical Stein extension operator as defined in Lemma 12.1 
then the native space of our approximant given by (|19|) is 



E:H T (J?;div) x H T -\f2) -> A/"*.(R d ) = tt T (R d : div) x iT"" 1 ^). 

Once again we can define interpolants with scaled kernels. In this case, we 
define the matrix- valued kernel 



&5 : = 



&T+1.5 

1pT-l,S 



(d+l)x(d+l) 



(24) 



where & T +i t s := (— Z\I-|-VV T )'0r+i,5 and the scaled basis functions are defined 
as in (0. Then the native space of the kernel $>s is given by 



A/a 



1 )=A^ T+M (M d )xA^ T ^ a (M d ), 
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with norm 



Ifullk 



(2tt 



-d/2 



llwlUVr+l,*^) Vr-l,5(w) 



do;. (25) 



We will need norm equivalence as stated in the following lemma. 
Lemma 3.1 For every S £ (0,<5 a ] where ip T +i and Vv-i generate H T+1 (M. d ) 
and H T ~ 1 (R d ) respectively, we have Af^ s (R d ) = 7V#(lR d ) and for every f e 
JV&(M. d ) there exist positive constants c$ and cq such that 

c 5 /2 H f Ha^*,, (R d ) < l|f||jV*(R<i) < c 6 II f II AT*,, («<*)• 

Proof For 8 < 1 we have 

(1 + 5*\\u\\iy > (5* + 5*\\u\\iy > g» (1 + . (26) 

With the definitions of the * and norms in (|2"2"|) and (|23|) . and with (fT5|) 
and (|16[) . we can see that 



H^pll^ T _i( Kd ) 



Ci 

Similarly, for S > 1 we have 



2 r2r+2| 



W*(R d )- 



l f ll^(R<2) ^ Jl fu H^ +1 

= ^"11 f II at*. (»<*)■ 



ll^pll^ T _l(» d ) 



For the lower bound if <5 < 1, we can simply use that (5 2 ||u>|| 2 < to see 

that 



lA/* 4 (R d ) ^ ^ 



< 



If 112 



* T+ i(R d ) 



II// 2 



— ^l|f|lA^*(K d )- 

If 6 > 1, with 

(1 + iMi 2 r = <^ 2t (<5 2 + <5 2 ii^ii2) r > s- 2t (i + 8 2 Mir, 

we reach 



1^ llX* s (m d ) - ^r^ 2r+2 



if 11 2 1 n f 11 2 



ll/pllAA Vvi (K d ) 



&e +2 iif" 2 



C-2 



Then the result follows with C5 
c 6 := (C 2 /C 1 ) 1 /2 max (M-+ 1 ). 



(Ci/C 2 ) 1 / 2 min(l,^+ 1 ) and 



10 



A. Chcrnih, Q. T. Le Gia 



Wc require one further result from [20j . 

Theorem 3.2 Let m G N and let tt C W 1 be a C m + 1 ' 1 smooth domain 
with outer normal vector n. For each f G H" l (i?) and g € H m+3 / 2 (9J?) wii/i 
Jg n g ■ ndS = 0, the nonhomogeneous Stokes problem HI)-© has a unique 
solution u G H" l+2 (l?) and p G iJ m+1 (J2) and 

l|u|| H >" + 2(f2) + |b||ff™+i(fi)/R < C (||f || H >"(J?) + l|g|lH'" + 3/2( 9i 7)) ( 27 ) 

Theorem 3.3 Lei r > 2+d/2. Assume that fl CM. d is a bounded, simply con- 
nected region with a CM' 1 boundary. Let f G H T " 2 (I?) and g G H T -V 2 (a/2) 
satisfy Jg^S ' ndS* = 0. Suppose the kernel 3? is chosen such that Af&(M. d ) = 
H T (K d ;div) x H T ' 1 (R d ). Then the approximation S x v given by (JTSJl to the 
Stokes problem (H])-© satisfies the error bound 

||v-S x v|| L2( „) < C^- 2 ||Ev-SxEv||^ (ffi<i) , (28) 

where h := max(/ii, /12) and the extension operator E is given by (|23p . 

Proof With the definition of the Sobolev space norms in (|5J) and assuming that 
we choose the representer for the pressure p such that i^i ( 1?) /r = |H|iji(r2) 
gives 

||v - S x v\\i, 2 (a) < ||u - S x u|j L2(fi) + ||p - SxP\\l 2 ({2) 

< ||u - s x u|j H 2 (fi ) + ||p - Sxp\\w(n) 

= ||u - S x u|j H 2(fi) + \\p - S x p\\m(n)/M 

< C\\Lv - LS x v|| L2(fi) + ||u - S x u\\ H3/ 2 {dn) , (29) 

where the last line follows from ([27)) applied to v — S x v with m = 0. We 
now extend the function v to Ev G H T (R d ) x H T - 1 (R d ) and note that the 
generalised interpolant S^v coincides with S^Ev. We now consider the two 
terms in the right hand side of (|2"9"]l separately. From (fit))), we have 

||Lv - LS x v|| L2(fi) < C^r'llEv - S x Ev||^ (Rd) . 

From (fTTj) , we have 

||u - S x u|| H 3/2 (af3) < Chl- 2 \\u - S x u\\ H , (n) . (30) 

Now we can write 

||u - S x u|| H -(fi) < ||u - Sxu||»r(fl) + Hp - S x p\\ H r-i {n) 

< ||E div u - S, Y Edivu|| aT(Rd;div) + \\Ep- S x Esp\\ H t- i(r") 

< C||Ev-S x Ev||^ (Rd) , 

and the stated result follows. 
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4 Multiscale symmetric collocation approximation 

We can now formally state our multiscale algorithm for the symmetric collo- 
cation solution of (HJ)- (J3J) which is stated as AlgorithmQ] To simplify notation, 
we write S^v = S^v and 3?i = &$ i and denote the mesh norms for the in- 
terior and boundary collocation points at level i as hi t % and /i2.i respectively. 



Algorithm 1: Multiscale symmetric collocation approximation to the 
Stokes problem 

Data: n: number of levels 



Xi : = {Xi i, X2,i}" =1 '■ the interior and boundary collocation 
points for each level i, with mesh norms at each level given by 
{hi^i, /i2,i}™— i satisfying c/ihi < h^i < fihi, where 
hi := max(/ii i, /i2,i) with fixed fj, G (0, 1), c G (0, 1] and 
hi sufficiently small 

{Si}"— l '■ the scale parameters to use at each level, satisfying 
Si = 3 /( T + 1 )^ p j s a g xe( j constant. 



Set M v = 0,f = f,g = g 
for i = 1, 2, . . . , n do 

With the scaled kernel solve the symmetric collocation linear system 



Result: Approximate solution at level n, M n v 
The error at level n, e n := v — M n v. 



We require a technical lemma regarding the error in the estimation of the 
velocity u. 

Lemma 4.1 Assume that u G H r (i?;div) with r > and let Edi v be defined 
by §§§ . Then we have the following bound 



begin 



{LS i v{x)) j = /i_ M (x), 1 < j < d, x e 
(S i v(x)) J . = 9i-ij(x), 1 < j < d, x G X 2 ,i- 



Update the solution and residual according to 



M;v = Mi_iv + SiV 
fi = fi_i - LSiV 

Si = Si~l - SiV 
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Proof With the definitions of the E^v^s and T operators, we have 





2 




E div u(u;) 




U) X 




2 du>- [ - 





< C 



E s Tu(w) 



- C\\ E sTu\\ L2 ( Rd) 

< C\\E s Tu\\ 2 mm 

< C\\Tu\\ 2 H1{n) 



< C\\u\\ 2 



where we have also used that the Es and T operators are bounded (Lemma 
l2"TD . 

The following theorem and corollary arc our main results on the conver- 
gence of the multiscale symmetric collocation algorithm for solving the Stokes 
problem. 

Theorem 4.2 Assume that fi and f , g satisfy the smoothness assumptions of 
Theorem \3.3\ Suppose the kernel $ is chosen such thatM^lW 1 ) = H r (IR d ; div)x 
H T ~ 1 (R d ) and define the scaled kernels by (|24p with scale factor Sj. Then for 
Algorithm]]] there exists a constant a\ such that 

l|Ee,-|| JV i J+1 (R<i) < aiHEej-illjv^fRd), 



(31) 



where ot\ is a constant independent of the point sets Xi, X2, ■ ■ ■ and Ee 3 is the 
extension operator for v defined in (|23[) applied to the error at level j defined 
in Algorithm]!] 

Proof With the notation Eej = Ev-MjEv = (u-Mj-EjiyU, p- MjE s p) T = 
(Edi v e u .j, Ese p j) T and with (|23|) . we have 



l Ee j'llA^* j+1 (R d ) - @i 



Edi V e U j (<■»■>) 



(1 + <5? + iIM|1) 



2 \T + 1 



=-h+h, 
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with 



h 



Ml2<; 



\>*\W>-, 



Edi v e u ,j(w) 



Edi v e u ,j(w) 



(1 + <£ +1 |M|2) T+1 + £se p ;» (l + S 2 3+l \Ml)" 



(1 + 5 2 J+l \\u\\i) T+1 + Ese'Au) (1 + 8 2 j+1 \\u\\lY 



da;, 



da;. 



For Ii, we can use that <5_j+ij|u;j|2 < 1, Lemmata 13.11 and 14.11 and Theorem 13.31 
to yield 

h < C (lle^-H^ + \\E s e pd \\l 2 
< C (||e UJ ||L 2 (fi) + ll e p,jlli 2 (j?) 

3 jllAf*(R d ) 



< C^ r - 4 ||Ee, 112 



< C 



6 2r+2 



|Ee 



-It -2 



Ee 



where we have also used that the interpolant is norm-minimal with respect to 
the 7V* j .(lR d )-norm. 

For I2, since 5j+i||li;||! — 1' we nav e 

(l + ( 5 2 +1 ||u,|| 2 ) r < 2^110,11)2^ < (25 2 +1 ||u,|| 2 ) T < 2> 2t (l + *?|M|I) T , 

and hence if /i, 6 < 1. we have 

h < CA i2r ~ 2 H Ee ill^ j (K t i) 
< C 2 /i 2r_2 ||Ee i _i|| 2 vr ^ (K£i) , 

where we have again used that the interpolant is norm-minimal with respect 
to the (K d ) -norm. The result follows with 

Corollary 4.3 There exist positive constants C3 and C4 such that 

\\v-M n v\\^ 2{ n)<C^(\\n\\^T {n) + \\p\\ H T- K a^ for n=l,2,... 

and 

||u - M„u|| L2(afi) < C 4 < (||u|| H T(fi) + \\p\\ht- for n = 1, 2, . . . 

T/ims i/ie multiscalc approximation M„v resulting from Algorithm]^ converges 
linearly to v in i/ie L^—norm in fl and on dfi if ai < 1. 
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Proof With Lemma 13.11 and Theorems 13.31 and 14.21 we have 



||v - M n v|| L2 (fi) < ||e„|| L2 (r2) 




< Ca n x (||u|| H -(«) + WpWht- i(n)) 



which proves the first result. For the second result, with ([3T)|) we can see that 



and the remainder of the proof is the same as for the first result. 



5 Condition number 

In this section, we present upper and lower bounds for the eigenvalues of the 
multiscale symmetric collocation algorithm for the Stokes problem. At each 
step of the multiscale algorithm, we need to solve a linear system resulting 
from the collocation conditions (f2"0f and ([2"Tj) on a set X = {xi, . . . ,xm}: 



Since the collocation matrix Ag is symmetric and positive definite, we know 
that the condition number is given by 



where X max (Ag) and A m i n (A,5) denote the maximum and minimum eigenvalues 



We will first need several technical lemmas concerning derivatives of the 
Wendland functions. 

Lemma 5.1 With spatial dimension d and smoothness parameter k = 2, 3, . . . 
let ipe^k be the original Wendland function. Then with x,y e R d andl < i,j < 
d and i ^ j , we have 



u - M„u|| L2(afi) < ||u - M„u|| H 3/2 (afi ) 

< Ch T 2 - 2 \\u-M n u\\ H r {n) 

< Chn~ ||Ee n ||jvi( K <j), 



A s b = (f g) 



k(A s ) 



^max(-A-t)) 



(32) 



of A s . 



<9yV^,fc( x _ y)l 



= 0. 
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Proof We recall that the Wendland functions are pieccwisc polynomials with 
support [0, 1] and we can write |21| 

2k+e 

Tpt,k(r)=J2 biri ' re IM ( 33 ) 

i=0 

and that the first k odd coefficients {&2i+i }i = o vanish. With the chain rule, 
where x - y = (x\ - yi, . . . , Xd — yd) and r = ||x- y|| 2 , we have 

W« - y) = ^-^-^ (v-2(r) - ^g(r) 

Using (j3"3")l . this last expression becomes 

/ w \ /2k+£ 2k+e 

^,fc(x-y) = ^ ^ EM*" 1)2 r'" 2 - £ ib, r*- 



\ i=2 i=l 
/2k+l \ 



E^ ri " 2 (34) 



i=l 



where 

_ r(c + n) 

[ )n - r(c) 

denotes the Pochhammcr symbol. Now the first three coefficients {bi}^ =1 arc 

bi = h= 

6a = (2 - 2)6 2 = 

63 - 0, 

since the first k odd coefficients of the Wendland polynomial are zero and 
k > 2. Hence we can write 

/2k+i \ 

%^,fe( x - y) = ( x i - - Vj) E ^ rl ~ 4 ' 



\ i=4 



and the result follows immediately. 



Lemma 5.2 FFif/i spatial dimension d and smoothness parameter k = 3,4, . . . 
lettjj£ t k be the original Wendland function. Then withx.,y £ M. d and 1 < 
d and i j, we have 

dijA 2 ip itk (x - y)| x =y = 0. 
Proof Once again employing the chain rule gives 

%^, fc (x-y) = ^ffi^ 
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With (|3"3")l we can rewrite this as 



2k+£ 2k+(. 

1—6 i—5 
2k+l 2k+i 2k+l 2k+? 



^Z\ 2 ^, fe (x - y) = {xi y -^p ^ 



2k+£ 

2Z~ b > 



r 'i-6 



i=4 i=3 i=2 i=l 

2k+i 



Since hi = C(i)bi, the first k odd coefficients {t^i+ilfLo are zero - Then we can 
determine other coefficients as 

b 2 = 30(6 2 - b 2 ) = 

6 4 = b 4 {60 - 15(3)2 + 12(2)3 - 7(1)4) = 

b 6 = b 6 (90 - 15(5)2 + 12(4)3 - 7(3) 4 + (1) 6 ) = 0. 

Hence since k = 3, 4, . . ., we can write 

2k+£ 



dijA 2 ifj eyk (x - y) = (xj - yi)( a; j - Vj) X! ^ 



i=8 



and the result follows immediately. 

Lemma 5.3 With spatial dimension d and smoothness parameter k = 2,3,... 
let ipg : k be the original Wendland function. Then with x, y £ K d and 1 < j < d 
we have 

%V^,fc( x -y)lx=y < °. 

and is independent of j. 

Proof With the chain rule, where once again x — y = (x\ — yi, . . . ,x<i — yd) 
and r = ||x — y||2, we have 

^, fc (x - y) = { ^-^ (Vij(r) - l4>)) + \4i(r). 

With Lemma [5TT1 the term in brackets is equal to zero when x = y. Using (|3"3"|) 
and noting that the first k odd coefficients are zero, this last term becomes 

2t+k 

i-2 



. T 
i=2 
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which means that the case of x = y, which is equivalent to r — 0, reduces down 
to 262 ■ Now combining positive terms into a generic constant C, we have from 
® 



b 2 = C(-l) 
= C 



k 

(-l) fe 



r(|-(fc-i)) 
= c(-i) k (-i) k - 1 < 0, 

where we have also used [9l 8.339.3] 

,'1 \ ^(-4)"n! 
r - - n = V^- 



2 J v (2n)! ' 

Lemma 5.4 With spatial dimension d and smoothness parameter k = 3, 4, . . . 
let i/je,k be the original Wendland function. Then with x, y S K d and 1 < j < d 
we have 

<%Z\ 2 ^,fc(x-y)| x=y < 0, 

and is independent of j. 

Proof With the chain rule, where once again x — y = [x\ — yi, . . . ,x<j — yd) 
and r = ||x — y 1 1 2 , we have 

x + Uflir) ~ ^JJ(r) + ^!(r) - ^g(r) + £$>(r) 

With Lemma 15. 1[ the first term in the previous expression is equal to zero 
when x = y. As before, we can write the second term as a series 

2k+e 

2 = 1 

Since k > 3, b ± = b 3 = b 5 = and equating coefficients gives 

b 2 = 6 2 (-6 + 3(l) 2 ) = 

b 4 = fc 4 (-12 + 3(3)2 - 3(2)3 + 2(1)4) = 0, 

which means we are left with 

2k+l 



Ev 1 



i=6 



18 



A. Chcrnih, Q. T. Le Gia 



Hence the case of x = y, which is equivalent to r = 0, reduces down to b 6 
which is given by 

b 6 = ((2) 5 + 2(3) 4 - 3(4) 3 + 3(5) a - 18)66 = 11526b. 

As before, combining positive terms into a generic constant C and noting that 
k = 3, 4, . . ., we have from [3] 



6 6 = C(-l)*(f) 
(-1)' 



C- 



r(i-(fc-3)) 
= c(-i) k {-i) k - 3 < 0. 

The next theorem gives a lower bound on the minimum eigenvalue of Aj . 

Theorem 5.5 Suppose the kernel 3? is defined by (|17[) and define the scaled 
kernel <f?s by (|24j) with a positive scaling factor 6. Then the smallest eigenvalue 
of the collocation matrix defined by (|20[) and (|21[) can be bounded by 



A min (A) > C 



qx 
8 



2t 



1x 



-d-2 



where the constant C is independent of the pointset X. 

Proof We follow the proof of [5J Theorem 4.1]. We will adopt the functional 
notation 



£u(v) 



(iv) l (x J ) for l^j^N, 1 < i < d, 

Vi(xj) for N + 1 < j < M, l<i<d, 



and 



Once again we will use the superscript y to denote that the functional acts 
with respect to its second argument. Then with f3 E M. dM , we need to show 
that 



d m 2 
E E &J&>,kSi*$, k *6(*-y)>c(yf) T q x d - 2 \\(3\\i (35) 

— 1 j,k— 1 



Now with the inverse Fourier transform, the left hand side of (|35[) becomes 
d M 



E E 0«A^eijSU**(x-y) 

— l j,k—l 



d M 

EE j 'i a '-j {u)) 

1=1 J=l 



$a(o;)da; 
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Now we define a second scaled kernel <& a by (J7J) with < a < 1 and a < 5. 
With and recalling that Vv-i satisfies (|T5)l gives 

= vv^i(M > ci, T _i (i + \\5u\\iy T+1 

-T + l 



CL.r-1 \- s 
>Cl,r-l(f) 9r " a (l + ||Hl2)" T+1 



— — - 1 — — ( T ) ^r-l,a(«). 
C2,t-1 

Since Vv+i satisfies f|15f) and with (JTSJ), we proceed similarly to get 
Ch>(w) = *7Z(5lj) = (\\duj\\ 2 2 I - S 2 lju: t ) ^(<5||u;j| 2 ) 



It 

2, ., ,T\ 



ci jT+ i ^-J (||aw|| 2 I - s 2 uu , 



2 



-J (||aw|||l - « 2 ww T ) (l + ||a^||i) 
> (^) 2T (||au,|||l - aW) vQ(a|M| 2 ) 



c 2,t+1 



* T+ i, a (w). 



C2,r + 1 

Since a/5 < 1, we have the following bound on 



and hence we have 

d M 



E E^A',fc^,**«(x-y)>c(^) E E A,ift', fe ^ )J -^, fe *a(x-y), 



and if we select a = < 1 such that we need only consider entries of the 
quadratic form corresponding to equal centres, with the definition of the scaled 
kernel in (0, this reduces to 

d M 

E E ft,j&^«^**( x -y) 

£,i'=l j,k=l 

>c(^) 2T ^ d E|f>L I" E ^^Vr + i(0)-a A _ 2 a i4 Vv-i(0) ] + 

f=i [j=i v j={i^}\i / 
E ft* ( E Qx%rtr-m 

j=N+l \ j={l:d}\% 
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since for interior centres wc have 



v 2 Ej=i:d\i d n A 2 i) T+1 {0) - d u ip T -i(0) for i = i\ 

M.jSi', it*!-" y)\j=k - 



z^Zy 2 Vv+i(0)-dii'VV-i(0) = for i^i', 

(36) 

with Lemmas 15.11 and 15.21 Similarly for the boundary centres 

£ f' y *(x - v)l ■ * = { ~ ^i=^V %V»r-i(0) for i = i', 

and then the result follows as 

d M 



E E %ft^4 fc *Kx-y)>^ q x d ~ 2 \W\\l 



i=l j,k=l 

with Lemmas 15.31 and 15.41 which give 



c := mm - E ^d^Vr+iCO) - <9iiVV-i(0), - E %^r-i(0) 

3={l:d}\i j = {l:d}\z 



> mm^\- E %4Vr+l(0)-fiKVT-l(0),- E Wr-l(0) 
j={l:d}\i j={-L:d}\i 



(d d 
- E d ]3 A 2 ^ T+1 (Q) - d n A-i(0), - E fyj VV-i(O) 

since Vv+i is a radial function and 9ji^> T _i(0) is independent of i from Lemma 
1531 

Our next result bounds the maximum eigenvalue A max (A,5). 

Theorem 5.6 Suppose the kernel &g is defined as in Theorem 15.51 Then if 
we assume that 

M < Chr d , (38) 

where M denotes the number of (interior and boundary) centres, then the 
largest eigenvalue of the collocation matrix constructed with &s defined by 
(|20p and (|21[) can be bounded by 

A max (A 5 ) < CS- d ~ 2 h- d , 

if 5 > 1 and by 

A max (A 5 ) < CS~ d ~ 6 h- d , 
if 5 < 1, where the constants C are independent of the pointset X. 
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Proof Using the notation from Theorem 15.51 together with Gershgorin's the- 
orem, we have 



d M 



|A m „(A,)-&j^*j(x,x)i<£; y, ie^ fc ^(x,y)i, \<i<d 

i'=i k=i 

which since $ is positive definite, using Lemmas 15.31 and 15.41 the defini- 
tion of the scaled kernels ([7]) and (f31)| and ([3T|) . if 8 > 1 

A max (A 5 ) < dM||^,.^.*ff(-,-)IU-(nxfl) 

(d d 
- ^ %^Vr+l,i(0) - 9llVr-l,5(0), - ^ d»VV-l,<5(0) 
i=2 j=2 

(d d 
- d n A 2 ip r+1 (0) - 9iiW-i(0), - XI %V>r-i(0) 
i=2 j=2 

where in the last step we have used that 

d n A 2 ^ T+hS (0) = r^4 2 rVr+i(0) < 5~ d ~ 2 d n A 2 ^ T+l (Q), 
since S > 1. If S < 1, we have 
Amax(Aj) < dM||^,.^.* 4 ( v )IU.(nxn) 

(d d 
- %^Vr+l(0) - 9uVr-l(0), - X 0J3VV-l(O) 

where in the last step we have used, for example, that 

dnW-M(O) = 5- d dn5- 2 4> T ~i(0) < r d - 6 9 u Vr-i(0), 
which completes the proof. 

We note that (|38[) will hold if, for example, the dataset is quasi-uniform, which 
means that hj/qj is bounded above by a constant. 

Now with ([52)1 and Theorems 15. 51 and 15.61 we obtain the following theorem 
where we write q 3 := qx d ■ 

Theorem 5.7 Suppose the kernel &s is defined as in Theorem \5.5l Then the 
condition number of the multiscale symmetric collocation matrix in Algorithm 
[7] is level- dependent and is bounded by 

, r \ 2r-d-2 



if 5 > 1 and by 



c [f.) h 



lT-d-1 



hj\ -- T ^ T (2r-d-6)-d-4 
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if 5 < 1. In the case of quasi-uniform datasets and hj < 1, these reduce to 

kj < ChJ 2r+2 . 

Proof The first two results follows with 5j = fii 1 - 3 /' T+1 ' and (|5"2"j) and Theo- 
rems [53] and [OjJ If the datasets are quasi- uniform, which means that hj/qj is 
bounded above by a constant, the final result follows by simplifying the first 
two expressions. 



6 Numerical experiments 

In this section, we present the results from applying the multiscale algorithm 
described in Algorithm [1] with Q = [0,1] 2 and v = 1 to the Stokes problem 
with exact solution given by 

, \ f 2cos(5xi)cos(2a;2)\ 

p{xi,x 2 ) = sin(3xi) sin(3x2) + C. 

This gives 

, \_( 58 cos(5xi) cos(2x2) + 3cos(3a;i) sin(3x2) \ 
t(xi,x 2 ) - ^ 145sill ( 5a;i ) sill (2 X2 ) +3sin(3.Ti)cos(3x 2 ) J 

and g equal to the restriction of u(x) to dfi . 

We use the C 8 Wendland radial basis function given by 

V> 6 ,4(||x||) = (1 - ||x||)+ (429||x|| 4 + 450||x|| 3 + 210||x|| 2 + 50||x|| + 5) , 

which is positive definite on R 2 and generates the Sobolev space if 5 - 5 (R 2 ) j21) . 
We use the same kernel for both ip T +i and VV— i- Consequently, in this case 
t = 4.5. Since d = 2, our approximate solution takes the form 

N I i/(922A0 T+1 (x - Xj) \ M /-5 22 V'r+l(x - Xj) 

Sxv(x) = atx,j -vdi 2 Aip T+ i (x - x.,-) + a ld d 12 i> T +i{x - xj) 

j=l \ -9iV't-i(x-Xj) / j=N+l \ 

n ( -vd 12 Aip T+1 (x - Xj)\ m / <9 12 7/>t+i(x - \ 

+ ^ Ct2,j ^5nZi?/v+l(x - Xj) + Q'2.j -9n?A r+ i(x - Xj) . 

i=l \ -^r-lCx-Xj) / j=N+l \ / 

We used five levels for the approximation, with iV equally spaced points for 
the interior point sets and 4(y/~N — 1) equally spaced boundary centres. The 
number of interior points, Nj, the number of boundary points, Mj — Nj, and 
the maximum mesh norms at each level, hj, are given in Table [TJ We note that 
the (maximum) mesh norms decrease by one half at each level and hence we 
select [i = \. For the scaling parameters, since r = 4.5, Algorithm Q] specifies 
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Level 


1 


2 


3 


4 


5 


N 


25 


81 


289 


1089 


4225 


M - N 


16 


32 


64 


128 


256 


h 


1/4 


1/8 


1/16 


1/32 


1/64 



Table 1: The number of interior and boundary points used at each level and 
the maximum mesh norm each level for the numerical experiment 



that 

with (3 constant. With the given value of hi in Table [IJ we select j3 such that 
Si = 10. This gives (3 = 18.779 and we use this to generate the other 5 values 
which are given along with the L2 and L x errors and condition numbers (re) 
collocation matrix in Table [5] The Li error was estimated using Gaussian 
quadrature with a 300 x 300 tensor product grid of Gauss-Legendre points 
and the error was estimated with the same tensor product grid. We used 
MATLAB for the calculations and worked with double precision. 



Level 


1 


2 


3 


4 


5 


«j 




10 


7.29 


5.33 


3.89 


2.84 


ll e u,i IIl 2 


(«) 


1.592e-02 


6.498e-04 


3.274e-05 


1.650e-06 


1.028c-07 


e ",i HLoc 


(fi) 


2.740c-02 


2.233c-03 


1.462e-04 


8.268e-06 


4.579c-07 


l|Ve pJ || L 


2 (fi) 


1.112e+00 


1.222c-01 


1.235c-02 


2.561c-03 


5.612c-04 


l|Ve p>J || ic 


o(«) 


4.209e+00 


3.338c-01 


1.048e-01 


3.650e-02 


1.211e-02 


Kj 


1.268e+09 


7.768c+ll 


2.573e+14 


6.879e+16 


1.674e+19 



Table 2: The scaling factors, approximation errors and condition numbers of 
the collocation matrices for the multiscalc symmetric collocation Stokes prob- 
lem example 
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